home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / gnu / glibc108.gz / glibc108 / glibc-1.08.1 / sysdeps / ieee754 / sqrt.c < prev    next >
C/C++ Source or Header  |  1992-03-12  |  3KB  |  121 lines
  1. /*
  2.  * Copyright (c) 1985 Regents of the University of California.
  3.  * All rights reserved.
  4.  *
  5.  * Redistribution and use in source and binary forms are permitted provided
  6.  * that: (1) source distributions retain this entire copyright notice and
  7.  * comment, and (2) distributions including binaries display the following
  8.  * acknowledgement:  ``This product includes software developed by the
  9.  * University of California, Berkeley and its contributors'' in the
  10.  * documentation or other materials provided with the distribution and in
  11.  * all advertising materials mentioning features or use of this software.
  12.  * Neither the name of the University nor the names of its contributors may
  13.  * be used to endorse or promote products derived from this software without
  14.  * specific prior written permission.
  15.  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR IMPLIED
  16.  * WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED WARRANTIES OF
  17.  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
  18.  */
  19.  
  20. #include <ansidecl.h>
  21. #include <errno.h>
  22. #include <math.h>
  23.  
  24. /* Return the square root of X.  */
  25. double
  26. DEFUN (sqrt, (x), double x)
  27. {
  28.   double q, s, b, r, t;
  29.   CONST double zero = 0.0;
  30.   int m, n, i;
  31.  
  32.   /* sqrt (NaN) is NaN; sqrt (+-0) is +-0.  */
  33.   if (__isnan (x) || x == zero)
  34.     return x;
  35.  
  36.   if (x < zero)
  37.     return zero / zero;
  38.  
  39.   /* sqrt (Inf) is Inf.  */
  40.   if (__isinf (x))
  41.     return x;
  42.  
  43.   /* Scale X to [1,4).  */
  44.   n = __logb (x);
  45.   x = __scalb (x, -n);
  46.   m = __logb (x);
  47.   if (m != 0)
  48.     /* Subnormal number.  */
  49.     x = __scalb (x, -m);
  50.  
  51.   m += n;
  52.   n = m / 2;
  53.  
  54.   if ((n + n) != m)
  55.     {
  56.       x *= 2;
  57.       --m;
  58.       n = m / 2;
  59.     }
  60.  
  61.   /* Generate sqrt (X) bit by bit (accumulating in Q).  */
  62.   q = 1.0;
  63.   s = 4.0;
  64.   x -= 1.0;
  65.   r = 1;
  66.   for (i = 1; i <= 51; i++)
  67.     {
  68.       t = s + 1;
  69.       x *= 4;
  70.       r /= 2;
  71.       if (t <= x)
  72.     {
  73.       s = t + t + 2, x -= t;
  74.       q += r;
  75.     }
  76.       else
  77.     s *= 2;
  78.     }
  79.  
  80.   /* Generate the last bit and determine the final rounding.  */
  81.   r /= 2;
  82.   x *= 4;
  83.   if (x == zero)
  84.     goto end;
  85.   (void) (100 + r);        /* Trigger inexact flag.  */
  86.   if (s < x)
  87.     {
  88.       q += r;
  89.       x -= s;
  90.       s += 2;
  91.       s *= 2;
  92.       x *= 4;
  93.       t = (x - s) - 5;
  94.       b = 1.0 + 3 * r / 4;
  95.       if (b == 1.0)
  96.     goto end;        /* B == 1: Round to zero.  */
  97.       b = 1.0 + r / 4;
  98.       if (b > 1.0)
  99.     t = 1;            /* B > 1: Round to +Inf.  */
  100.       if (t >= 0)
  101.     q += r;
  102.     }                /* Else round to nearest.  */
  103.   else
  104.     {
  105.       s *= 2;
  106.       x *= 4;
  107.       t = (x - s) - 1;
  108.       b = 1.0 + 3 * r / 4;
  109.       if (b == 1.0)
  110.     goto end;
  111.       b = 1.0 + r / 4;
  112.       if (b > 1.0)
  113.     t = 1;
  114.       if (t >= 0)
  115.     q += r;
  116.     }
  117.  
  118. end:
  119.   return __scalb (q, n);
  120. }
  121.